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We introduce fractal liquids by generalizing classical liquids of integer dimensions d — 1,2,3 
to a non-integer dimension dj. The particles composing the liquid are fractal objects and their 
configuration space is also fractal, with the same dimension. Realizations of our generic model 
system include microphase separated binary liquids in porous media, and highly branched liquid 
droplets confined to a fractal polymer backbone in a gel. Here we study the thermodynamics and 
pair correlations of fractal liquids by computer simulation and semi-analytical statistical mechanics. 
Our results are based on a model where fractal hard spheres move on a near-critical percolating 
lattice cluster. The predictions of the fractal Percus-Yevick liquid integral equation compare well 
with our simulation results. 


The liquid state, an intermediate between gas and 
solid, exhibits short-ranged particle pair correlations in 
isotropic shells around a tagged particle [l|. While the 
shell structure is lost for gases, solids exhibit long-ranged 
correlations and anisotropy. Particle correlations are ac¬ 
cessible by experiments (1[| and contain valuable infor¬ 
mation about the particle interactions. A fundamental 
task in theory and computer simulation of the liquid state 
is to predict particle correlations for given interactions. 
In this respect, one particularly successful approach is 
liquid integral equation theory [l], . 


Molecular and colloidal liquids can be restricted to one 
or two spatial dimensions (5lJ[6| by confining them on 
substrates 0, at interfaces , between plates 0, 
in channeled matrices |ll| or optical landscapes 
The thermodynamic properties change accordingly [If 
and the systems are called one-dimensional (ID) or two- 
dimensional (2D) liquids. While these classical exam¬ 
ples of confinement involve integer-dimensional configu¬ 
ration spaces, there are also cases where liquids are con¬ 
fined in porous media [l3 - [l7| . or along quenched poly¬ 
mer coils [l8| . and where the configuration space ex¬ 
hibits non-integer (fractal) dimension at suitable length 
scales. Much work has been devoted to understand¬ 
ing the limiting cases of low and high densiW namely 
the motion of a single particle on a fractal |l9l42.3l | and 
the structure of a fractal aggregate itself, corresponding 
to an arrested high-density particulate system [iiSi. 
In all these situations, particles interact in the embed¬ 
ding integer-dimensional space and interaction energies 
depend on Euclidean particle distance. 


In this letter, we break new ground by considering 
fractal particles in a fractal configuration space, both 
of the same non-integer dimension. An alternative 
model in which the particle dimension differs from the 
configuration-space dimension is briefly mentioned near 
the end of the letter. Figure [T] is a snapshot from one 



Fig. 1: (Color online) Snapshot of a Monte Carlo simulation 
where 300 fractal particles are located on a percolating cluster 
shown in blue. The dimension of each particle and of the 
cluster is di = 1.67659, and the chemical distance diameter 
of each particle is cr = 300a (c.f.. Fig. mi. Red, black and 
gray pixels are occupied by particles. Red: Particle center 
regions (chemical distance to a particle center is less than 
O.lcr). Black: Particle rim regions (chemical distance to a 
particle center is more than 0.4a and less than 0.5cr). Gray: 
Regions between particle centers and rims. Every pixel in 
the figure corresponds to one vertex in the simulation. One 
quarter of the whole simulation box is shown, containing here 
72 particle centers. 


of our Monte Carlo (MC) simulations. The interaction 
between two particles is not described by a pair poten¬ 
tial in Euclidean space, but instead inherits fractal char¬ 
acter from the particles and their configuration space. 
Unifying the two complementary fields of classical liquid 
state theory and fractal structures, we refer to the ex- 
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plored, dense disordered phases as ’fractal liquids’. As 
opposed to continuum-mechanical descriptions of fluid 
flow on fractals or liquids forming fractal stream 

networks here we calculate the correlations of in¬ 

dividual fractal particulates that constitute a liquid. We 
develop basic concepts for statistical mechanics of frac¬ 
tal liquids by presenting MC simulations along with an 
accurate semi-analytical fractal liquid integral equation 
approach. Our model system is a fractal hard sphere liq¬ 
uid, the analogue of the simplest generic model liquid in 
integer dimensions. In the present proof of principle, we 
restrict our studies to the liquid state. Determination 
of the full phase diagram of fractal particles in fractal 
space is a task that should be tackled in future work. In 
particular, it will be interesting to study whether or not 
symmetry breaking and crystallization occur in fractal 
spaces. 

Fractal liquids are not only interesting from a theoret¬ 
ical perspective, but are also found in nature. Porous 
media filled with phase-separated binary liquids such as 
water and oil appear in natural oil or gas reservoirs, and 
are produced during hydraulic fracturing in oil drilling. 
If the oil droplets exceed the porosity length scale in size, 
they adapt to the fractal shape of the void space. Such 
droplets constitute fractal ’particles’ in a fractal config¬ 
uration space set by the porous medium. While the in¬ 
dividual water and oil molecules are confined to the void 
space, it is important to note that the mesoscopic wa¬ 
ter or oil droplets are not necessarily confined: If the 
size of the void space greatly exceeds the droplet size, 
then the droplets move in a (quasi) unbounded configu¬ 
ration space of their own non-integer dimension. Here we 
study the collective bulk-phase thermodynamics of such 
droplets, which are the fractal liquid constituents. Other 
possible model systems include different mesoscopic ’par¬ 
ticles’, such as floppy slime molds [sj or droplets con¬ 
fined to fractal polymer aggregates in a gel. 

To model a fractal liquid we first need a fractal¬ 
dimensional configuration space with a suitable distance 
measure. A viable choice is a large set of vertices on 
a discrete lattice for a fractal lattice liquid simulation. 
From a regular square lattice with a number W of ver¬ 
tices and periodic boundary conditions in embedding 2D 
Euclidean space, we randomly remove vertices until V 
vertices are left. According to percolation theory 0,[33, 
in the limit IT oo there is a critical percolation thresh¬ 
old Pc = T/IT at which a fractal-dimensional percolating 
cluster extends to infinite length scales. For the square 
lattice, a threshold of pc = 0.5927... has been reported 


Throughout this work, all distances I between vertices 
are evaluated as shortest distances on the percolating 
cluster, using our own implementation of what is essen¬ 
tially equal to Dijkstra’s algorithm [s^. The distance 
between two vertices is I = a ii and only if these vertices 
are connected by a bond in one of the two Cartesian di- 



Fig. 2: (Color online) Lattice schematic. From an 8 x 8 square 
lattice with periodic boundaries, 17 vertices have been re¬ 
moved. Circles are the remaining vertices and black lines are 
the bonds of length a. The chemical distance between the 
orange and the blue filled circle is I = 8a along any of the 
three paths indicated by red, blue and orange arrows. The 
corresponding Euclidean distance is r = 2a. 


rections of the embedding space. Illustrated in Fig. [5J 
this distance measure is commonly referred to as ’chemi¬ 
cal distance’ and the corresponding non Euclidean metric 
is known as the ’taxicab metric’. 

The lattice in Fig. [5] features the essential properties 
of the lattices in our simulations, except of their much 
larger size of IF = 3000 x 3000, F = 5 355 000. The 
ratio VjW = 0.595 is 0.4% larger than pc, which ensures 
the presence of a percolating cluster in the finite system. 
We retain only those vertices that belong to the percolat¬ 
ing cluster. All other vertices, either in smaller clusters 
or isolated, are deleted. An average count of 3 x 10® con¬ 
nected vertices remains. As shown in Fig. [31 the average 
number P(/) of vertices with chemical distance I from a 
tagged vertex grows as a non-integer power of 1. 

The so-called ’spreading dimension’ di = 1.67659 ... of 
the percolating cluster is the relevant dimension here. It 
has been determined by MC simulations which agree 
with our result P{1) oc Note that di differs from 

the more commonly reported fractal cluster dimension 
df = 91/48 [20|, [ 34 ], [37| , because chemical distance scales 
as a non-integer power of Euclidean distance. 

We define a fractal particle as the set of vertices that 
are within a distance Z < cr/2 from a center vertex. There 
is exactly one center vertex per particle and we refer to cr 
as the particle diameter. Every particle occupies at least 
one vertex pair with distance Z = cr at all times. Particles 
obey a no-overlap constraint, prohibiting configurations 
in which any two particle centers have a distance of Z < cr. 
We choose cr = 300a, with the consequence that the par¬ 
ticles themselves are fractal objects (c.f. Fig. [3]). As 
the snapshot in Fig. [T] shows, the particles are highly 
anisotropic in the embedding 2D space. However, in 
fractal space and taxicab metric, the particles and their 
no-overlap interactions are perfectly isotropic. In spite 
of their polymorphy in embedding space, all particles 
are indistinguishable for the purpose of our simulation. 
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Fig. 3: (Color online) Fractal clustersize scaling. Solid curves 
indicate average numbers, P{1), of vertices with chemical 
distance I from a tagged vertex, in the largest cluster re¬ 
maining from a lattice of initially W = 3000 x 3000 ver¬ 
tices. Results are for various numbers IF — F of removed 
vertices, with ratios VjW as indicated. Upper dashed line: 
P(Z) = 4 Z/a, corresponding to a 2D lattice. Lower dotted 
line: P(Z) = 1.51 (Z/a)°'®^ corresponding to a 1.67659- 
dimensional lattice. Our simulations are for VjW = 0.595. 


Hence the particles are fractal-dimensional analogues of 
monodisperse hard spheres, and one may refer to the sim¬ 
ulated system as the 1.67659-dimensional hard sphere liq¬ 
uid. 

We probe the thermodynamic equilibrium state by MC 
simulation. Particles are moved one at a time, and ev¬ 
ery move proceeds as follows: After deleting a random 
particle, we pick a random vertex vi globally from the 
cluster, which is a candidate to become a vertex at the 
rim of the displaced particle. We then pick a random 
vertex V 2 at distance I = a from vi , and a random vertex 
V 3 at distance Z = tT/2 from both vi and V 2 - Vertex is 
the candidate to become the displaced particle’s center 
vertex. If there is a center vertex of another particle at 
distance I < a from vs, then the move is rejected and 
the original particle restored. Otherwise, the move is ac¬ 
cepted. 

The simulation starts with a random configuration of 
potentially overlapping particles with diameter a = 10a. 
In the initial simulation stage a is gradually inflated to 
300a, which facilitates finding a dense configuration with¬ 
out overlaps. After the initial inflation and equilibration 
the production stage is entered, and the average num¬ 
ber N{1) of particles with center-to-center distance I is 
recorded. Binning the function N{1) with a bin width of 
5a reduces scatter in the data. Inside a fractal control 
volume, the fractal packing fraction (j) is measured as the 
average number ratio of vertices occupied by particles to 
the total number of contained vertices. The control vol¬ 
ume is a set of vertices with distance I < 1500a from a 
center vertex. It is randomly moved during the simula¬ 
tion and contains at least one vertex pair with distance 
I = 3000a at all times. We compute the chemical distance 



Fig. 4: (Color online) Equations of state for hard spheres in 
dimensions from cZ = 1 to d = 3. Thick blue curve: Exact 
(Tonks gas / PY) result for d = 1. Thick red curve: Kolafa- 
Rottner (3^ virial series with 15 coefficients for d — 2. Thick 
green curve: Carnahan-Starling equation of state for d = 3. 
Thin black curves are the PY virial pressure results for di¬ 
mensions as indicated. 

distribution function 

g{l)=c X N{1) X (1) 

with the Z-independent factor c chosen such that g{l) 1 
for large Z. Function g{l) is the analogue of the radial dis¬ 
tribution function g{r) of isotropic liquids in integer di¬ 
mensions [J. Our simulations require long runtimes due 
to the numerically expensive distances calculations. The 
fractal cluster’s ramified structure severely complicates 
any performance improvement, except for parallel exe¬ 
cution of statistically independent runs and subsequent 
averaging. 

The Ornstein-Zernike equation for a homogeneous and 
isotropic liquid in integer-dimensional Euclidean space 
reads 

g(r) = 1 + c(r) + n J d‘^r'[g{r') - 1] c{\r - r'\), (2) 

where n is the particle number density, c(r) is the di¬ 
rect correlation function [l|, r is a d-dimensional vector 
with Euclidean norm r = |r|, and d‘^r' is an infinitesimal 
volume element at position r'. In conjunction with the 
no overlap constraint g{r < a) =0 and the approxima¬ 
tion c(r > ct) = 0, Eq. ([2]) constitutes the Percus-Yevick 
(PY) integral equation for cZ-dimensional hard spheres 
[39|, which can be systematically derived by functional 
Taylor expansion [l|. 

We solve the PY equation by means of a spectral 
solver (^ . Our numerically efficient algorithm for the 
convolution-type Eq. dZD is based on the Hankel trans¬ 
form pair 

Nd/2 

•^( 9 ) = \d/2-i J (3) 

0 

i-d/2 y 

/W = TTTm / dq (4) 

( QiTt ) J 

^ ^ 0 
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for a d-diniensional isotropic function /, sampled on a 
logarithmic grid (dll-ld^. In Eqs. m and 0, Jdl2-l{x) 
is the Bessel function of the first kind and order c?/2 — 1. 
Since Jd/ 2 -ii^) is analytic with respect to both x and 
d, the solution can be carried out formally also for non¬ 
integer dimensions. Resulting pair-correlations and ther¬ 
modynamic properties represent the analytic continua¬ 
tions of standard PY theory with respect to the dimen¬ 
sion. 

The capability of liquid integral equations to pre¬ 
dict thermodynamic properties in arbitrary dimension 
is exhibited in Fig. 01 featuring equations of state 
for hard spheres in various dimensions from d = 1 
to d = 3. Thin black curves in Fig. |4] represent 
the d-dimensional hard sphere reduced virial pressure 
p/nksT = 1 -f 4)g(r = a~^) with Boltzmann con¬ 

stant ks, absolute temperature T, packing fraction (j) = 
and g(r) calculated in the PY 
scheme. Thick curves in Fig. 0] are exact or nearly ex¬ 
act reference solutions for d = 1,2 and 3. For d = 1, 
the PY result reduces to the exact (Tonks gas) solu¬ 
tion p/nksT = 1/(1 — (j)). The good agreement of the 
PY-scheme with the numerically accurate result from 
Ref. [ 3 ^ for d = 2 verifies the fidelity of the PY scheme 
for 1 < d < 2. Comparison to the Carnahan-Starling 
equation of state p/nksT ss (1 -|- (/ -|- (/^ — (j)^)/{l — (j))^ 
for d = 3 reveals a decreasing PY-scheme accuracy for 
higher dimensions at large values of (j). 

Figure 0] shows the functions g{l) extracted from our 
simulations with 300 and 150 particles. Also shown are 
the PY-scheme solutions g{r) for d = 1, d = di = 1.67659 
and d = 2 and for (j) = 0.487 and 0.266. With the values 
of (j) determined in the simulation, the PY scheme is free 
of adjustable parameters. The PY-scheme solutions for 
d = di = 1.67659 agree best with the simulation data, 
which is most clearly seen for the higher density liquid 
(Fig. 01 upper panel). Note that the chemical distance 
range in Fig.0]is the same as in Fig. 01 where the fractal¬ 
dimensional cluster scaling is validated. 

We have also simulated 2D hard disks with center-of- 
mass positions confined to a fractal configuration space 
of dimension d < 2 (results not shown). This constitutes 
an alternative fractal liquid model system, where parti¬ 
cle dimension differs from configuration space dimension. 
Pair correlations in such liquids differ distinctly from the 
predictions of the fractal-dimensional PY scheme. This 
can be rationalized by noting that the PY-scheme equa¬ 
tions are based on one and only one measure of distance, 
which is incapable of describing both a fractal-space cor¬ 
relation and an embedding-space interaction. 

In conclusion, we have simulated a fractal liquid in 
thermodynamic equilibrium and demonstrated that the 
measured particle correlations are well predicted by a 
fractal liquid integral equation. Our approach may be 
easily applied to other configuration spaces with differ¬ 
ent fractal dimension. Particle interactions beyond the 
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Fig. 5: (Color online) Symbols: Chemical distance distribu¬ 
tion functions g{l) from our MC simulations. Curves: PY 
solutions g{r), for dimensions and (fractal) packing fractions 
as indicated. Top and bottom panels are for 300 and 150 
simulated particles, respectively. 


simple no overlap constraint should soon be studied, in¬ 
cluding short-ranged attraction which should result in 
liquid-vapor demixing, and also long-ranged repulsion. 
We have reported here on a monodisperse model system, 
but interaction-polydispersity is already implemented in 
the liquid integral equation solver for arbitrary dimen¬ 
sion. 

Fractal liquids lead the way for a plethora of new 
fundamental research. At present it is unclear which 
types of thermodynamic phases exist in fractal dimen¬ 
sions and phase diagrams await to be outlined. Field 
theories like density functional theory could be general¬ 
ized to non-integer dimension. Mode coupling theory for 
the kinetic glass transition of three-dimensional liquids 
in porous media (3 and in arbitrary dimensions 
relies on static structure input, which can be calculated 
using fractal liquid integral equations. Future studies 
should cover time-resolved dynamics of fractal liquids in 
and out of equilibrium. Transport coefficients of fractal 
liquids may be studied, requiring an account for fractal 
hydrodynamics. A promising application for fractal liq¬ 
uid theory is the prediction of thermodynamic properties 
of microphase separated liquids in porous media as en¬ 
countered in natural oil and gas reservoirs. 
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